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Recently, we developed and implemented the bond propagation algorithm for calculating the par- 
tition function and correlation functions of random bond Ising models in two dimensions.— The 
algorithm is the fastest available for calculating these quantities near the percolation threshold. In 
this paper, we show how to extend the bond propagation algorithm to directly calculate thermody- 
namic functions by applying the algorithm to derivatives of the partition function, and we derive 
explicit expressions for this transformation. We also discuss variations of the original bond propa- 
gation procedure within the larger context of YVY-reducibility and discuss the relation of this class 
of algorithm to other algorithms developed for Ising systems. We conclude with a discussion on the 
outlook for applying similar algorithms to other models. 



I. INTRODUCTION 

For nearly 80 years Ising models have given valuable 
insight into phase transitions and critical phenomena 
in magnets, alloys, and many other systems. Random- 
bond Ising models (RBIMs) in particular are often used 
to study frustration and spin-glass behavior, and they 
are closely related to neural networks and information 
theory— In Ref. Q], we presented a numerically exact^I 
algorithm for computing the partition function and cor- 
relation functions in a class of 2D Ising models, which 
works for any planar network of Ising spins with arbi- 
trary bond strengths but without applied fields in the 
bulk. Applications include random-bond Ising models 
(RBIM) (including ±K disorder, Gaussian disorder, site 
dilution, and bond dilution) and geometric frustration as 
in the case of triangular Ising antiferromagnets. Here, we 
show how to extend the algorithm to directly calculate 
thermodynamic quantities such as the internal energy. 

Since its introduction in 1925, many exact results have 
been found for regular 2D Ising models (i.e., those with 
translational symmetry in two directions). There has 
been continued interest in random-bond Ising models, 
which have more complexity because of the need to av- 
erage not only over all spin configurations, but also over 
all configurations of bond strengths. Approximate ana- 
lytic methods can give misleading or conflicting results, 
so numerical calculations play an important role in the 
field. 

Numerical methods for 2D RBIMs have developed 
along several lines. Onsager's original solution using 
operator techniques generalizes to the fermion-network 
method of Merz and Chalker— . The Ising partition func- 
tion is related to the number of dimer coverings (per- 
fect matchings) of a modified graph, and hence to the 
problem of computing a determinant or Pfaffian^; this 
approach was applied to RBIMs by Saul and Kardar^i 
and by Galluccio et a&. The fastest of these algorithms 8 



takes of order 0(N 3 ^ 2 ) time for a network of N spins. 
Algorithms in this class have the disadvantage that one 
must first map to fcrmions or dimcrs before the model 
can be solved. 

A different line of development involves the Y-V and 
V-Y transformations for the partition functions of Ising 
modcls^^ operating directly in the spin basis, without 
the need to map to fermion or dimer models. In a Y-V (or 
star-triangle) transformation, three Ising spins connected 
by bonds to a central spin can be converted into the same 
three Ising spins but without the central spin, now con- 
nected by mutual bonds, in such a way as to preserve 
the total partition function of the system. The reverse 
V-Y transformation can also be done so as to preserve 
the partition function. For networks which are YVY- 
reduciblc, one can then compute the partition function 
Z at a given temperature using a suitable sequence of 
such transformations. Colbourn et alr^ suggested us- 
ing the Feo-Provan-^ YVY reduction method for gen- 
eral networks, which takes 0{N 2 ) time where N is the 
number of nodes in the network (i.e., spins). Frank and 
Lobb invented the bond propagation algorithmic (a form 
of YVY reduction) for 2D resistor networks which takes 
0{N Z / 2 ) time for square lattices (grid graphs)^ They 
also suggested an extension to Ising models. Recently, we 
developed and implemented the Ising bond propagation 
algorithm^ The algorithm executes in 0(N 3 ^ 2 ) time for 
most planar networks of interest, and in NhiN time for 
dilute models near percolation^ making it the fastest 
method for computing the Ising partition function and 
correlation functions near the percolation threshold. In 
this paper, we derive the transformations necessary to 
implement the algorithm on derivatives of the partition 
function, allowing for fast, direct, and exact calculation 
of thermodynamic quantities. We also discuss the rela- 
tion of this class of algorithm to graph theory as well as 
to other methods for Ising systems. 

The outline of this paper is as follows. In Sec. |TT] we 
introduce the models to which the algorithm is applica- 
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ble. In Sec. IIIII wc derive the network transformations, 
most notably the Y-V and V-Y transformations, that are 
the building blocks of the bond-propagation algorithm 
(BPA) for both resistor networks (i.e., Gaussian mod- 
els) and the partition function of Ising models, including 
new forms of the transformations which have improved 
numerical stability in the latter case. In Sec. IIVI we de- 
scribe some variations of the algorithm. In Sec. [V] wc 
show that the BPA can be used to compute derivatives of 
the Ising partition function directly, and we derive these 
transformations for the first derivative with respect to 
temperature. In Sec. I VII we discuss practical issues such 
as the strengths and limitations of the BPA, as well as 
the outlook for application of similar algorithms to other 
models. 



II. MODELS 

YVY reduction is applicable to both Gaussian models 
(such as resistor networks) and Ising models. We define 
a 'Gaussian model' in statistical mechanics as one whose 
action S({4>}) is bilinear in the fields 0, which may be 
real, complex or Grassmannian, 

e^W^expV^.^.. (i) 

ij 

This action also applies to dynamical models that are di- 
agonal in frequency space, such as electrical LCR net- 
works involving inductors, capacitors and/or resistors 
where the represent admittances (complex conduc- 
tances), and non- interacting bosons, fermions, phonons 
or magnons on lattices with hopping amplitudes K+j. 

An Ising model action is a quadratic form in spin vari- 
ables <7j = ±1 with a symmetric kernel Kij, 

e s(M) = e -/m(M) = exp Y^K^i (2) 

(ij) 

where (ij) indicates a sum over pairs of spins. K = f3K 
are dimensionless Ising couplings, K are the physical 
couplings with dimensions of energy, and (5 = 1/T is 
the inverse temperature. This includes various types of 
random-bond Ising models (RBIMs) such as those with 
'binary' couplings (ify = ±Kq), Gaussian couplings, and 
zero couplings (dilute Ising models). Such models are of- 
ten used to study spin glasses and dilute magnets. 

In the Ising system, our goal is to calculate the parti- 
tion function, 

z = £y<w), ( 3 ) 

M 

correlation functions such as 

(a i a j ) = Z- 1 J2a i a j e s ^\ (4) 



and/or thermodynamic functions such as 



U = -Z-i§ = Z-^nW})e S({ ° }) - (5) 

For Gaussian models such as resistor networks, we aim 
instead to calculate the effective resistance Rij between 
two points in the network. 

Computing the partition function of a many-body sys- 
tem is usually a very difficult problem that involves sum- 
ming over a number of configurations that is exponential 
in system size. Even if this is re-cast into the form of 
a high-temperature or low-temperature series, it is still 
usually the case that one needs to sum over a large num- 
ber of graphs, or count a large number of graph embed- 
dings, using some amount of brute force. Treatments 
based on successive elimination of spins are often limited 
because integrating out spin degrees of freedom tends to 
generate longer range couplings or multispin interactions 
such as 3-spin couplings. Another approach is real-space 
decimation, which is approximate because one has to dis- 
card most interactions beyond a certain level of complex- 
ity. 

There are some nice exceptions: for Gaussian models 
and Ising models on planar networks^ it is possible to 
reduce an entire network by successive Y-V and V-Y 
transformations, local moves which never induce multi- 
spin interactions. Thus throughout the whole procedure 
the system can be described just by pairwisc interactions; 
therefore, it is useful to adopt the language of graph the- 
ory and represent it by a network, with nodes (vertices), 
edges (bonds), and edge weights (couplings). The next 
section describes transformations that can be used in the 
reduction of such a 'statistical mechanics network'. 



III. STAR-MESH AND MESH-STAR 
TRANSFORMATIONS 

A. Gaussian models 

For pedagogical reasons, wc first review Gaussian mod- 
els (which include resistor networks as a special case). We 
show how to eliminate any node of a Gaussian model via 
a star-mesh transformation regardless of its coordination 
number. Consider the network in Fig. [T] described by 
an action S((j)o, 4>\, . . . , 4>n) containing bilinear couplings 
Kij(/)*(j)j as well as a constant 'free energy' term F. The 
variable 4>o can be eliminated by Gaussian integration, 
leading to a new effective action S'((f>i, . . . , 4>n) with pa- 
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FIG. 1: Elimination of an iV-degree node, illustrated for 
N = 3. Diagonal couplings K\\, K22, K33 correspond to loops 
(edges connecting a node to itself) and have been omitted. 



rameters F' and K' 



{d(j> ) expS(<f>o,<pi,..-,<j>N) 
f (# ) exp U + Kijtffa ] 

( -1 N N ( 



-Koo 



(N N 
t=i j=i 

exp S'(<f>i, . . . ,(/>n) 

rf(3fty)d(S<ft) 



(6) 



and £ = 1 in the case where 



where J (d<p) = J 
4> are complex variables; similar expressions exist for real 
4> and Grassmannian <p. The integration generates addi- 
tive changes to the free energy, SF = F' — F, and to the 
couplings, SK i:j = K' H - K id : 



SF = Cm 



-1 



K 00 ' 



on 



i,j = 1,2,3,... ,N. 



(7) 
(8) 



Eq. ([5]) contains only algebraic operations — addition, 
multiplication and division — and it is homogeneous, i.e., 
invariant under multiplication of all if's by the same con- 
stant. In fact, Eq. ([8]) corresponds to row and column 
subtractions on the the matrix Kij ; the elimination of 4>q 
by Gaussian integration is closely related to the proce- 
dure of Gauss elimination in linear algebra. 

The above derivation is valid even if the variables 
. . . 4>n are coupled to external fields or additional vari- 
ables; these extra terms simply cancel out on both sides 
of Eq. ©. 

Node elimination corresponds to the "star-mesh trans- 
formations" illustrated in Fig. [2] and defined as follows. 
An N -degree star-mesh transformation eliminates a node 
of coordination number N, introducing N{N— 1)/2 edges 



N 



Star 



Mesh 
(empty) 



K 01 
• O 



K01 



a-;., 




Transformation name 

Elimination of 
free spin 

Elimination of 
dead end, pendant, 
or dangling bond 

Series reduction, 

decoration-iteration 

transformation 



Y-Delta, Y-nabla, 
wye-delta, 
star-triangle 
transformation 



FIG. 2: Star-mesh transformations for various N. Circles 
represent nodes and lines represent edges; the filled node is 
the one being eliminated. 



between its neighbors. A mesh-star transformation is the 
inverse, which does not always exist. 

For Gaussian models the star-mesh formula is 



-K t0 K Qj /Ko 



,3 = 1,2,3,. 



(9) 



In the special case of a resistor network, the off-diagonal 
elements of the kernel are equal to the conductances, 
= Gij, and the diagonal elements are given by a 
sum rule K vl = -^jLi^iji 

becomes G'^ = G i0 G j / 'X)j=i G oj 
responds to the familiar series reduction of conductances 
For N 



so the star-mesh formula 
For N = 2, this cor- 



3, it reduces to the V-V transformation, 



G' 12 = GoiG 02 /(G 01 + G02 + G 03 ) (cycl.), (10) 

where 'cycl.' indicates additional equations in which the 
indices 1,2,3 are cyclically permuted. Eqs. (flU)) can be 
inverted to give the V-Y transformation, which is most 
easily written in terms of the resistances R = 1/G: 

R01 = JM 2 /(^i + R i2 + ^23) (cycl.)- (11) 

The equations for node- voltage analysis of a planar elec- 
trical network are the same as those for loop-current anal- 
ysis of the dual network, which is obtained by interchang- 
ing node voltages Vi with loop currents J M , and conduc- 
tances Gij joining adjacent nodes with resistances 
separating adjacent loops. Thus the similarity between 
Ecis. (fTU)) and (fTTj) is no coincidence. 

For computer implementation it may be preferable or 
necessary to supplement the above formulas by formulas 
for special casesi^, involving zero couplings (opens) or 
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infinite couplings (shorts). Assuming that floating-point 
overflows never occur, the following rules are sufficient: 

YDelta(oo, Gi, G 2 ) = (0, G 2 , Gi) (12) 
DeltaY(0, Gi, G 2 ) = (oo, G 2 , Gi) (13) 
DeltaY(oo, G X ,G 2 ) = (Gi + G 2 , oo, oo). (14) 

The Gij can be complex numbers representing admit- 
tances at a particular frequency, in which case the trans- 
formations involve complex arithmetic. 



B. Ising models 

For Ising networks, one may attempt to derive a star- 
mesh transformation in a similar way, integrating out a 
spin Co to generate an effective action S' with parameters 
F' and K' to be determined (see Fig. EJ: 

e %,» I ,-.w) =e s'K™>w) \/<ri...a N 

tr =±l 

(15) 

<T = ±1 

(16) 



Because oi, . . . , ctjv can each take the values ±1, there 
are 2 N equations that need to be satisfied. Due to spin- 
flip symmetry, only 2 N ~ 1 of these are independent. For 
N > 4 the system of equations is over determined, with 
no solution. For example, for N = 4, there are 8 indepen- 
dent equations in 7 unknowns F' , K' 12 , A{ 3 , A{ 4 , K 23 , 
K' 3i , K' 3i ; the true effective action S' contains a 4-spin 
interaction of the form K' 123i aia 2 a 3 ai that does not fit 
within the network formalism. It is futile to keep track 
of such multispin terms because their number can grow 
to 0{2 N ). 

For N = 3, however, there are 4 independent equations 
in 4 unknowns F' , K[ 2 , K' 13 , A 23 , so a Y-V transforma- 
tion exists^ such that 



E 



7 F-\-Kai<T0<7 1 + ^02t J "o£72 4--ft'o3 (T CT 3 



= 2e F cosh(A i<7i + K 02 a 2 + K 03 <r 3 ) 

_ e F' +K' 23 a 2 <J3+K'. 31 (j 3 ai+K' 12 (j 1 (72 (17) 

Writing SF = F' — F and the new couplings K'^ in terms 
of the old couplings Kio, 

2 cosh (Koiat + K 02 a 2 + K a3 a 3 ) 

= e SF + K '-23 a -i<y3+K! il a 3 a 1 +K' 12 <y 1 a 2 \ja\ (T 2 (7 3 

'oo = 2 C Osh( + A 01 + K 02 + K 03 ) = e SF+K^+K 31 +K' l2 

oi = 2cosh(-A i + K 02 + A 03 ) = e SF+K 23 -K 31 -K' 12 
a 2 = 2cosh(+A 01 - A 02 + A 03 ) = e'^W-^ 
a 3 = 2cosh(+A 01 + A 02 - A 03 ) = e 5F - K 23- K ' 3 i+ K 'i2 

'SF = lna Q 1 / i a 1 1/i a 2 1 / i a 3 1/ \ 



K' 23 = lna 1/4 ai 1/4 c 



■ *a 2 1/4 a 3 1/4 (cycl.) 



(18) 



where ai are functions of Aoi as above. 

For computer implementation it is preferable to write 
the Ising action as 



const -z JJ% (1 "°' <0 '' 



-<^)/2 



(19) 



where 2 = e F and fcjj = e 2Kij . In this representation 
the Y-V transformation is 

YDelta(fc 1 , fc 2 , fc 3 ) = (fc 23 , fc 31 , fc 12 ; <5z) where 
5z = 1 + kik 2 k 3 
Z\ = k\ + k 2 k 3 (cycl.) 

b = ziz 2 z 3 /Sz 
k 23 = b/zi (cycl.) (20) 

where we have written fcoi = and <5z = e SF (and omit- 
ted primes on k 23 etc.). The transcendental functions 
(cxp, log, cosh) have been replaced by algebraic func- 
tions (+, x, ^J~), which arc quicker to evaluate and 
have fewer complications arising from multivalucdncss. 
Furthermore, infinite ferromagnetic couplings K = 00 
('shorts' in the language of Ref. [HI ) are readily repre- 
sented by k = 0. 

The V-V" transformation is equivalent to a Y-V trans- 
formation in the dual representation using Syozi's 'cyclic 
change of lattices—: 

DeltaY(fc 2 3,fc3i,fci 2 ) = (ki,k 2 ,k 3 ) where 
Pi = dual(/c 23 ) (cycl), 
(P23iP3l>Pl2) = YDelta(pi,p 2 ,p 3 ), 

fei = dual(p 23 ) (cycl.). (21) 

where the duality relation is the Mobius transformation 
dual(fc) = fcp However, this approach does not give an 
expression for the free energy change 8 A,— and we have 
found that a direct implementation of Eq. (|2ip is sus- 
ceptible to roundoff error. The scheme below has better 
numerical performance: 



x\ = 2fc 3 ifci 2 (l - k 2 



(cycl.) 



2/1 = 1 + k 31 2 k 12 2 - k 12 2 k 23 2 - k 23 2 k 31 2 (cycl.) 



ki = xi/(yi + v) or (y t - u)/xi (cycl.) 

5z = 1 + fcifc 2 A: 3 (22) 

The two expressions for fci are formally equivalent; the 
choice depends on which one is numerically more stable 
(i.e., avoids subtraction of similar quantities). 

Eos. ([20]) and {22]) arc preferable to those in Refs. [Ill 
and [l|, because they each contain only one square root, 
and have been optimized to reduce roundoff error. Both 
transformations are l-to-2 mappings by virtue of the 
square roots, whose signs can be chosen arbitrarily^ For 
the V-Y transformation this reflects the Z 2 symmetry of 
the action: for a given Y, one can obtain an equivalent 
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Y by flipping the signs of Co, -Koij -K02, and -K03 simulta- 
neously. There is no clear physical interpretation of the 
doublc-valucdness of the Y-V transformation. 

Setting one or more of the fcj to 1 in Eqs. [20l gives the 
star-mesh formulas for N = 2,1, 0: 

^ = 2: fci 2 = ± fc f 2 , ^ = l + feife (23) 

iV = 1 : 2 = 1 + *i (24) 
AT = : z = 2 (25) 

The N = case (integrating out a free spin) simply in- 
creases the constant 'free energy' term by In 2, i.e., by 
one bit. 

The special cases below are useful for computer im- 
plementation. They are equivalent to Eqs. (fT4|) . with G 
replaced by K: 

YDelta(0,fci,fc 2 ) = {1MM) 
DeltaY(l,fci,fc 2 ) = (0,/c 2 ,fci) 
DeltaY(0,fci,fc 2 ) = (fcifc 2 ,0,0). (26) 

IV. THE BOND PROPAGATION PROCEDURE 

The transformations derived in the previous section 
can be used to reduce a network (i.e., eliminate all or 
nearly all its nodes) and hence to calculate the partition 
function or effective coupling(s). Although an inefficient 
method, gaussian models can be reduced simply by elim- 
inating the nodes one by one, in any order, using the N- 
star-mesh transformation, Eq. (|9|). In contrast, for Ising 
models a general star-mesh transformation does not ex- 
ist. However, certain networks can be reduced using only 
-/V = 3,2,1,0 star-mesh and mcsh-star transformations; 
these are known as Y VY-reducible networks. The bond- 
propagation algorithm is an efficient example of YVY- 
reduction. 



A. The original bond-propagation algorithm 

Gaussian models and Ising models on square lattices^ 
can be reduced by the bond propagation procedurei^ 
depicted in Fig. [3J A single bond propagation move is 
illustrated in Fig. [3]^a). First, a V-Y transformation is 
performed on the upper left V of Fig. E^a). This intro- 
duces one new node. The new node is now effectively 
shifted so as to replace the node to its lower right, a 
"move" which docs not change the topology of the net- 
work. Then, a Y-V transformation is used to convert 
the lower right Y into a V, removing a node. In this 
way, any diagonal bond can be "propagated" into a di- 
agonally adjacent plaquette. The prescription originally 
implemented in Ref. [j~3| is shown in Fig. EJb). Starting 
from the upper left corner in Fig. [3Jb) , a series reduc- 
tion is used to convert the corner into a diagonal bond. 




(a)A single bond propagation move 




(b)Latticc reduction 



FIG. 3: The bond-propagation algorithm in the original form 
invented by Frank and Lobb. (a) A single bond propagation 
move consists of a Y-V transformation followed by a V-Y 
transformation, (b) To reduce a finite square lattice, first a 
corner node is eliminated by series reduction, producing the 
first diagonal bond, which is then propagated diagonally down 
and to the right until it is "absorbed" at the opposite edge. 
Other corner nodes are eliminated in the same way. 



Then using successive bond propagation moves, this di- 
agonal bond can be moved diagonally down and to the 
right until it annihilates at an edge with open boundary 
conditions. 

For resistor networks, repeated bond propagation 
moves reduce the network to a single string of conductors 
in series, which is easily reducible to one effective con- 
ductance. For an Ising model, the effective coupling is 
related to the correlation function between corner spins: 
((Jii&nn) = tanhA'njvjv- One may also wish to calcu- 
late the partition function in the Ising model; in order to 
do this, one must collect the contributions to the parti- 
tion function during every transformation of the network. 
At the end of the calculation, after the last spin or node is 
eliminated, one has the free energy; the penultimate step 
gives the effective Ising coupling or effective resistance. 

The number of operations necessary to accomplish this 
is of the order of 0(L 3 ). Ref. HH showed empirically that 
by taking advantage of the early termination of bond 
propagation on dilute lattices, it is possible to achieve 
0(L 2 \nL) computational time scaling near the percola- 
tion threshold. 



B. Useful variants 

Figure [4] illustrates a variant of the algorithm for a 
(finite) triangular lattice. Here, by embedding the trian- 
gular lattice in a square lattice, we see that one third of 
the bonds can be considered "diagonal" to begin with. If 
these diagonal bonds are only propagated out when it is 
necessary to make space for other diagonal bonds, then 
this procedure reduces the size of the lattice by one unit 
in every direction. Alternatively, the diagonal bonds can 
all be propagated out, and then the BPA can proceed as 
with any square lattice. 

Any planar network can be embedded in a square or 



6 




FIG. 4: Reduction of a triangular lattice of side L to one of 
side L — 1 using ^L 2 bond propagations and L series reduc- 
tions. 




FIG. 5: A cylinder, when viewed in perspective from one end, 
is equivalent to an annulus, which is planar and hence YVY- 
reducible. 

triangular lattice (e.g., by inserting 'opens' and 'shorts' 
which do not change the partition function), to which 
the BPA can then be applied. A network with cylin- 
drical boundary conditions is equivalent to an annular 
network, which is planar (see Fig. [5]). Networks with 
toroidal boundary conditions, however, are non-planar. 

Self-averaging quantities can be calculated efficiently in 
a strip geometry. For the rectangular strip in Fig. [6j use 
a series reduction on the upper left corner to produce the 
first diagonal bond, which is propagated down and to the 
right until it is removed from the system. This creates 
two new corners in the upper left, to which the same 
procedure can be applied. As corners are eliminated, 
this leaves a string of bonds with couplings Kf, iff, and 
so on. For sufficiently large n, the couplings K^,K^ +1 , 
etc., come from the same distribution. These couplings 
may be used to determine the spin correlation length £. 
In this way, the BPA is easily adapted for the calculation 
of self-averaging quantities in strip geometries (Fig. [6j. 
For a rectangular Lx M strip, its time requirement scales 
as 0(LM 2 ). 

Instead of using the BPA YVY reduction procedure, 
efficient parallel algorithms may be developed by direct 
YVY reduction on every plaquctte simultaneously. The 
algorithm shown in Fig. [7] represents the same arithmetic 
operations performed in a different order. Starting with 



FIG. 6: Bond propagation for strips. From the asymptotic 
distribution of K„ one can infer quantities such as the asymp- 
totic effective resistance per unit length as a function of strip 
width M (for a random resistor network), or the correlation 
length £m for a random-bond Ising model of width M such 
that (eri,ier£,i) ~ e~ i/?JU for large L. 




3x3 triangular lattice Honeycomb 2x2 triangular 



FIG. 7: Highly parallel algorithm for reduction of a triangular 
lattice. Starting with a triangular lattice of A sites, perform 
V-Y transformations to insert B sites. On the resulting hon- 
eycomb lattice, perform Y-V transformations to eliminate all 
the A sites. This gives a triangular lattice that is one unit 
smaller in every direction. 




duality duality 




infinite couplings 
along boundary 



FIG. 8: A variant of Fig. [7] using duality and Y-V transfor- 
mations only. Thick lines represent infinite-strength bonds. 



a triangular lattice, one performs a V-Y transformation 
on every upward-pointing V, thus producing a honey- 
comb lattice. One then performs a Y-V transformation 
on every downward-pointing Y, producing a triangular 
lattice one unit smaller than the original. 

Using the 'cyclic change of lattices' introduced by 
Syozii^, it is even possible to devise an algorithm involv- 
ing only Y-V and duality transformations, as in Fig. [8] 
The dual of a planar network is derived by drawing a 
dual bond cutting across each bond of the original net- 
work, with a bond strength given by the duality relation 
dual(G) = (for resistors) or dual(fc) = (for Ising 
bonds). 

The latter two algorithms for YVY reduction are very 
suitable for parallelization. In Figs. [Jj and [U the 0(L 3 ) 
operations can be distributed among L 2 processors so 
that the entire network can be reduced in only O(L) time. 



C. Reduction of networks with external terminals 

When the BPA is used to compute the effective cou- 
plings between an even number of 'external' or 'terminal' 
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1 

o 



3 

6 

(a) 





(b) (c) 
FIG. 9: Irreducible 2-, 3-, and 4-terminal networks. 



spins (2,4,6,...), one must take care not to eliminate 
these terminal spins. Given any two-terminal planar net- 
work, it is always possible to integrate out all the 'inter- 
nal' spins, reducing the network to a single bond between 
the two terminals (Fig. 9(a)). However, for a general 4- 



tcrminal planar network, it is not possible to reduce the 



network beyond the stage illustrated in Fig. 9(c) which 
has one internal spin. It may be interesting to find irre- 
ducible networks with larger numbers of terminals (sec, 
e.g., Ref.[Ii). A rectangular L x M strip (L > 2M) with 
terminals all along its short edges can be reduced to a 
2M x M strip by successive applications of the proce- 
dure in Fig. [TU1 this may be thought of as a decimation 
procedure which preserves the form of the Ising Hamil- 
tonian. 

For a triangular network with one terminal at each 
corner, the algorithm in Fig. [7] leads to an irreducible 
three-terminal network (Fig. 9(b) I, from which the ef- 
fective two-terminal couplings K^i, K\2 may easily 
be found. For uncorrelated random-bond models, this 
approach produces effective couplings for three disorder 
realizations for the price of one. 



ical differentiation. 

The quantity U = — ■^klnZ(P) can be computed from 

-fiZU = ^S({a})e 5({CT}) 
M 

= Yl ( ]C K v a i a i) ex P Y KijVtO- 3 (27) 
L nodes 



M 
nodes 
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V. BOND PROPAGATION FOR CALCULATING 
DERIVATIVES OF THE PARTITION FUNCTION 

The algorithm described above (and its variants) al- 
lows us to calculate the numerical value of the partition 
function Z(/3) at one value of /3. To compute other ther- 
modynamic quantities such as the mean energy U, one 
approach is to calculate Z at closely spaced values of 
(3 and perform numerical differentiation. This is not 
appealing because while the algorithm produces near- 
machine-precision values for Z, this would yield only 
approximate values for U. We now show that it is in 
fact possible to extend the bond-propagation algorithm 
to compute derivatives of Z directly, i.e., without numer- 



FIG. 10: 'Contraction' of an L x M rectangular strip pre- 
serving nodes along short edges, where L > 2M. Dotted lines 
denote bonds that are removed by Y-V transformations; bold 
lines indicated bonds that are propagated out. Each curved 
arrow indicates a sequence of bond propagations. 



In order to use the BPA to calculate this quantity di- 
rectly, we develop Y-V and V-Y transformations which 
preserve this quantity. First, consider a 'Y' with terminal 
spins 01, 02, 03 and central spin a®. In fact, we can solve 
a more general expression, in which the bond strengths 
in the prefactor can differ from those in the exponent, 
such that the 'Y' is described by 8 parameters A, L i, 
L02, L 3 and F, K m , K 02l K 03 : 



C= J2 E + A + L ^ a « + L " 2CT2(T ° + L 03Wo) e s ^ +F+Ko ^ +Ko ^ +K ^° . (28) 

[fl,2,3,...] "O 

In the above expression, the quadratic terms in the prefactor and exponent not involving ctq have been collected in 
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the functions c ex t(ci, (72, ... ) and S^xtO^i, (72, ■ ■ ■ )■ Factorizing and evaluating the sum over Co, 



[ CT 1,2,3 



= ^ e Sc X t+Fl 2(c cxt + ^)cOsh(if lCri+A"o2Cr2 + A"o3Cr3) 
r l ^ v " 



[ CT 1,2,3 



1st term 

2( L 01 a! + L Q2 a 2 + L 03 a 3 ) sinh (Koxa! + K 2O-2 + K 03 a 3 ) }. (29) 



2nd term 

We wish to find the new parameters A' L23, L% 1: L f 12 and i 7 ", ^12 °f the equivalent V such that 

C= ( CcXt +A ' + L 23^ + L 3lWl + L' 12 CT l( 72) e S - +F ' +if -^2 ff 3+X3 1 -3^+K{2- 1 -2 



["'1,2,3,. 



gSoxt+Fj ^ c ^ + j^ e &F+K' 23l 72<J3+K' 31 <j 3l 7 1 +K' 12 <7 1 <j. 2 

1 '- S ' 

.2, 3, ..-J n ^ 1 

' ' ' 1st term 

+ (<L4 + £ 23 (7 2 (73 + ^31(73(7! + L^G^e^+^^+i^i^r+K^^ | 



(30) 



2nd term 



where 6 A = A' — A and SF = F' — F. Comparing Eqs. (J^HJ) and (|3T))) suggests equating the two terms inside the 
braces separately. Equating the first term leads to 

2cosh (Kaxai + K 02 a 2 + K 03 a 3 ) = e * F + K 'v'*"*+K'»i"*'i+Kl*<n"» Vai, a 2 , <t 3) (31) 

which is identical to Eq. (jT5J) . Equating the second terms in the braces in Eqs. I|29p and ([3D)) gives 

2{L 01 a 1 + L Q2 a 2 + L 03 a 3 ) sinh (K 01 a 1 + K Q2 <r 2 + K 03 a 3 ) 

= (SA + L' 23 a 2 <j 3 + L' 31 a 3 a 1 + L[ 2 a 1 a 2 )e SF+K -^ +K ^^ +K '^^ V*i , a 2 , a 3 . (32) 
We can now use Eq. (fT8]) to rewrite Eq. ([32]) as 



2(L iCi + L 02 a 2 + L 03 a 3 ) sinh (ifoioi + K 02 a 2 + K 03 a 3 ) 

= (SA + L' 23 a 2 a 3 + L' 31 a 3 a x + L' 12 crio- 2 )2cosh [K Q1 ai + K 02 a 2 + K 03 a 3 ) V<7i, a 2 , a 3 (33) 
(SA + L' 23 a 2 a 3 + L 31 o 3 o\ + L' 12 aia 2 ) = (ioi<7i + L 02 a 2 + L 03 cr 3 ) tanh (Koxcri + K a2 a 2 + I<o 3 a 3 ) (34) 

f SA + L' 23 + L' 31 + L' r2 = u = (+L 01 + L Q2 + L 03 ) tanh(+/v 01 + K Q2 + K 03 ) 

[SA + L' 23 - L' 31 - L' 12 = ui = {-Loi + L 02 + L 03 ) tanh(-Jf i + K 02 + K Q3 ) (cycl.) 

jSA = (uo + ui +u 2 + u 3 )/4 

\ L 23 = ( u o + ui - it 2 - u 3 )/4 (cycl.) 

where u% are functions of Loi and K$i as defined above. 



I 

This Y-V transformation can be inverted to give a V- where 
Y transformation. First, compute Kqi, Kq 2 , Kq 3 using 

Eq.©. Then let fc = coth(A' 01 + K 02 + K 03 ) 

(c - Cl +c 2 +c 3 )L' 2;i +cyci. 1 c i = coth(-A 01 + J\ 02 + A' 03 ) (cycl.). 



^ -C0+C1+C2+C3 ' 

>i = ^-^ + c 3) l 23+ ^ (cycLfeqs. ©,©,09,© can be used in the bond- 

(37) propagation algorithm (or any other YVY reduction 
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scheme) to compute the mean energy U of an Ising net- 
work directly, without numerical differentiation. Note 
that since the transformations which preserve Z are not 
the same as those which preserve its derivative Z 1 ', obtain- 
ing U = —Z'/Z requires two separate bond propagation 
calculations: one to calculate Z; the other to calculate 
Z' . It may be possible to calculate the heat capacity and 
other higher derivatives of Z in a similar manner. 

VI. DISCUSSION 

We have found that bond propagation generally gives 
answers to high accuracy (close to machine precision), 
provided that one uses a careful implementation that 
avoids overflow and underflow and minimizes roundoff 
error, especially in the case of large systems, low temper- 
atures, or frustration. We discuss specific cases below. 

A. Ferromagnetic models 

The algorithm works rather straightforwardly for ferro- 
magnetic Ising models where all the are positive. The 
Y-V and V-Y transformations should be implemented 
in a way that minimizes roundoff error, as in Eqs. ([20]) 
and (|22|) . For large systems in the ferromagnetic phase 
[L > 10 3 ), the partition function Z may cause floating- 
point overflow, which can easily be dealt with by storing 
F = In Z instead. The couplings themselves kij may also 
underflow, in which case it may be necessary to use an 
alternative representation to store the couplings and per- 
form the transformations, or to change the order of the 
bond propagations. 

B. Dilute models 

Applying bond propagation to dilute models contain- 
ing zero couplings may generate infinite couplings, as 
pointed out by Frank-Lobb. We have confirmed that this 
can be successfully dealt with using Eqs. (fT4|) or (|26| . 

C. Frustration 

All couplings Kij in a physical Ising model must be 
real- valued. The Y-V transformation, Eq. (f!2t)]) . preserves 
this property, but the V-Y transformation, Eq. ([^)) . can 
generate complex-valued couplings. This happens if the 
argument of the square root is negative (v 2 = yi 2 — x\ 2 < 
0), i.e., 

(1 + k 23 k 3 i + k 31 k 12 + k 12 k 23 ) 
x (1 + k 23 k 31 - hik 12 - k 12 k 23 ) 
x (1 - fc 2 3&3i + k 31 k u - k 12 k 23 ) 
x (f - k 23 k 31 - k 31 k 12 + k 12 k 23 ) < 0. (39) 



This "frustration inequality" can be taken as the defi- 
nition of a "frustrated V"^ Frustration occurs in V's 
where all 3 couplings have similar magnitudes and 1 or 
3 or them are antiferromagnetic (Kij < 0, kij > 1). The 
BPA still works for frustrated systems provided that com- 
plex arithmetic is used. Despite the occurrence of com- 
plex intermediate couplings, we have found that the value 
of the partition function emerging from the entire calcu- 
lation is real, as expected, with a small imaginary part 
arising from roundoff error. 

The models typically studied in the literature are 
random-bond Ising models with bond strengths picked 
independently from the same distribution, P(K). The 
two most common models are the 'binary' (±Ko) RBIM, 
with P(Kij ) = (1 -p)S(Kij - K ) +pS(K ij + K ) where p 
is the concentration of antiferromagnetic couplings, and 

the Gaussian RBIM, with P(Kij) oc exp-± . 

The BPA works well for the Gaussian RBIM. For the 
binary RBIM, certain disorder realizations contain per- 
fectly frustrated plaqucttes; these cause the BPA to gen- 
erate infinite and zero couplings, which then lead to in- 
determinate results. This is due to singularities in the 
function represented by the bond propagation procedure 
(the composition of all the Y-V and V-Y transforma- 
tions). These can be avoided using perturbations away 
from perfect frustration, including roundoff error itself. 
Such perturbations may be relevant in a renormalization- 
group sense, changing the universality class; if the scaling 
flow is weak, it may still be possible to obtain accurate 
results this way. Introducing perturbations this way is 
acceptable as long as the systematic error it produces 
does not exceed the random statistical error from dis- 
order averaging, which is an inevitable part of numerical 
calculations on RBIMs. At very low temperatures round- 
off error may become an issue for both types of RBIM. 



D. Comparisons with other algorithms 

Gaussian models (including resistor networks) possess 
./V- star- mesh transformations for all N > 1, so there are 
a large number of algorithms for reducing them. Below 
is a comparison of the number of multiplications ( x ) and 
divisions (-=-) required to perform various algorithms on 
Gaussian models (e.g., resistor networks) on an L x L 
lattice: 

• Cramer's rule: 0(L 2 l) (x), 0(1) (-=-) 

• Full Gauss elimination: 0(L 6 ) (x) , 0(L 2 ) (-=-) 

• Transfer matrices!^: 0(L 4 ) (x) , 0(L 3 ) (-=-) 

• Banded Gauss elim.: 0(L 4 ) (x) , 0(L 2 ) (-=-) 

• Nested dissection^: 0(L 3 ) (x) , 0(L 2 ) (-=-) 

• Bond propagation! 3 -: 0(L 3 ) (x) , 0(L 3 ) (-=-) 
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The last two methods are the fastest, both taking 0(L 3 ) 
time overall. The nested dissection algorithm devel- 
oped by George^ and by Lipton et aL— is a hierarchical 
"divide-and-conquer" method that exploits prior knowl- 
edge about the connectivity of the network. It is quite 
a general idea and it can be used with doubly periodic 
boundary conditions, as well as for lattices in dimensions 
d > 2, for which the time requirement is 0(L 3d ~ 3 ). Its 
drawback is the amount of bookkeeping required. On the 
other hand, bond propagation, while slightly more costly 
in terms of divisions, is simpler to implement and espe- 
cially to parallelize. (See Sec. lIVBl ) For dilute networks 
near the percolation threshold, bond propagation is even 
faster, taking 0(L 2 ha.L) time, because many propagat- 
ing bonds terminate early. In the context of Gaussian 
models, bond propagation can be viewed as an efficient 
linear algebra method for reducing certain pcntadiagonal 
matrices while preserving their sparsity. 

For Ising models, polynomial-time methods are only 
known to exist in two dimensions. Below are operation 
counts including addition (+), subtraction (— ), multipli- 
cation (x ), division (-=-), and square roots (^/) for selected 
methods: 

• Transfer matrices in the spin basis^: 0(2 L ) op- 
erations (+x) and memory; can be generalized to 
higher-dimensional lattices 

• Fermion network method 3 -: 0{L A ) operations (+x) 

• Pfaffian elimination with nested dissection^: 0(L 3 ) 
operations (+ x -=-) 

• YVY reduction by Feo-Provan algorithm^: 0(L 4 ) 
operations (H — x -j- for general graphs 

• YVY reduction by bond propagation^: 0(L 3 ) op- 
erations (H — x -=- 

Due to the lack of a general TV-star-mesh transformation 
for Ising models, nested dissection cannot be applied di- 
rectly to the Ising network. Rather, one must first map 
the Ising problem to a dimcr problem on the Kasteleyn 
'terminal lattice' which has 4 nodes for each Ising spin. 
In two dimensions, the dimer problem maps to a Pfaf- 
fian elimination problem, to which nested dissection can 
then be applied. In contrast, the bond propagation al- 
gorithm can be performed directly in the Ising network 
representation, with no need to map to dimers or fermion 
models. 



E. Limitations 

YVY reduction techniques (of which the bond- 
propagation algorithm is one) for statistical mechanics 
models rely on two criteria: (a) YVY-reducibility of the 
graph on which the model is defined and (b) the existence 
of Y-V and V-Y transformations that preserve some as- 
pect of the model at hand. 



We first discuss criterion (a). The Robertson-Seymour 
theory of forbidden minors and obstruction sets (see 
Ref. [2CJ and references therein) implies the following re- 
lationships among various sets of graphs: 

{Forests} 

C {Outcrplanar graphs} 

C {Series-parallel-reducible graphs} 

C {Planar graphs} 

C {YVY-reducible graphs} 

C {Linklessly cmbeddable graphs}. (40) 

In particular, this shows that all planar graphs are YVY- 
reducible. However, it is also true that nearly all 'suf- 
ficiently large' non-planar graphs with crossing number 
2 or greater are not YVY-reducible. This means that 
YVY approaches such as bond propagation are unable 
to completely reduce 3D lattices or even 2D 'pyrochlorc' 
lattices (lattices of edge- or corner-sharing tetrahedra). 
For example, a large enough 3D network with simple cu- 
bic symmetry has as a minor the Petersen graph KG^^- 
Since KG$ t 2 is a forbidden minor for YVY reduction, 
the simple cubic lattice is not amenable to the BPA, nor 
to any algorithm based purely on YVY reduction. Since 
general 3D lattices contain the simple cubic lattice as 
a minor, they are not YVY-reducible either. As a fur- 
ther example, although diagonal bonds can be made to 
propagate on, e.g., a tiling of 3-space by truncated oc- 
tahedra, this 3D lattice is not reducible by bond propa- 
gation. Of course, even for networks that are not com- 
pletely YVY-reducible, it may be useful to perform par- 
tial YVY-reduction before resorting to more expensive 
computational techniques. 

We now turn to criterion (b). In Secs. lIIII andlVlwe de- 
rived Y-V and V-Y transformations for Gaussian models 
and Ising models. These models have the special prop- 
erty that 3-spin interactions are forbidden by symmetry. 
For most other models in statistical mechanics, elimi- 
nating a 3-coordinatcd spin generates 3-spin interactions 
which cannot be described within the network formalism. 
For example, consider the Ising model in the presence of 
magnetic fields hi (which may be uniform or random), 
described by the action 

e s(M) = cxp ( Kijcno-j + hm) • (41) 

(ij) » 

Let us integrate out the center spin of a Y and find the ef- 
fective action for the resultant V by equating Boltzmann 
weights for each of the 8 configurations of the outer spins 
((Ti, 0-2,0-3 = ±1): 

^""^ e S(<TQ,<Tl, <T 2 ,g 3 ) _ g S'(ffl,ff2,<f3) ^42) 
CO 

If we restrict ourselves to the form of Eq. (|4"Tj) , there are 
only 7 possible terms: S' = F + h\a\ + h^oi + /13173 + 
K230-20-3 + K31O3O1 + Ki20~i<J2- Therefore S' must con- 
tain a 3-spin interaction .K'1230'10'20'3, and it cannot be 
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written in the same analytic form as the original action: 
a Y-V transformation does not exist. 

Table U shows a 'feasibility study' for Y — V transfor- 
mations for the g-state Potts model, with couplings Kij 
that tend to align adjacent spins in the same 'direction' 
in spin space, 

e S(M) = cxp KyS^ , ^£{1,2,..., q}Vi. (43) 

(«) 

a generalized Potts model (where every bond is described 
by a matrix Py representing the weight factors for every 
combination of directions of two adjacent spins), 

e s « ff » =exp^f^, Ji6{l,2,.„,{}V«. (44) 

m 

and the g-state clock model (where q > 3), 

e S(M) = CX p^ ATy cos ^(a t - a,). (45) 

<«> 

For these models there are, again, more equations than 
parameters, so exact Y-V transformations do not exist. 

In the limit q — ► oo the clock model becomes the XY 
model. This suggests that Y-V transformations do not 
exist for XY models, nor for other models with continu- 
ous spins (e.g., classical Heisenberg models). 

This paper has focused on statistical mechanics, but 
YVY-reduction and the BPA also have important ap- 
plications in combinatorics and operations research. Ex- 
act Y-V and V-Y transformations exist for the shortest 
path problem, for the max flow (or min cut) problems^ 
for counting spanning trees and perfect matchings,— , 
and even for 2D foams in mechanical equilibrium^!. 
Such transformations do not exist for network reliability 
problcm a 2233 , nor for equilibrium flow problems with ar- 
bitrary potential-flow relations^ (i.e., electrical circuits 
composed of elements with arbitrary non-linear current- 
voltage characteristics). Thus, the former group of prob- 
lems can be solved on planar networks in polynomial 



time, whereas the latter group cannot. In Table HT1 all of 
the above problems are classified according to the maxi- 
mum degree of star-mesh transformations that they sup- 
port. 

The above examples have been for systems defined 
on undirected graphs, i.e., 'symmetric' bonds. YVY- 
reduction (and the BPA) may also be applicable to di- 
rected graph problems such as dynamical Ising cellu- 
lar automata.— Note that in systems where exact Y-V 
transformations do not exist, approximate Y-V transfor- 
mations may still have some use.— 



VII. CONCLUSION 

In conclusion, we have shown that the bond propaga- 
tion algorithm (BPA) may be applied not only to the 
partition function and correlation functions in the Ising 
model, but that it can also be applied to thermodynamic 
quantities via transformations which preserve the first 
derivative of the partition function. Similar transforma- 
tions may exist for higher derivatives as well. Since the 
BPA is a form of YVY reduction, it is limited to net- 
works which arc YVY-reducible. This includes planar 
graphs, and the BPA may be applied to Ising models on 
planar networks as well as planar resistor networks. In 
general, 3D lattices are not YVY-reducible, and so the 
BPA is not amenable to 3D resistor networks or Ising 
models in 3D. 



VIII. ACKNOWLEDGMENTS 

It is a pleasure to thank R. Fisch, E. H. Goins, 
C. J. Lobb, and F. Merz for helpful discussions. This 
work was supported by Purdue University and Research 
Corporation (YLL) and by the Purdue Research Foun- 
dation (EWC). EWC is a Cottrell Scholar of Research 
Corporation. 



1 Y. L. Loh and E. W. Carlson, Phys. Rev. Lett, 97, 227205 
(2006). 

2 H. Nishimori, Statistical Physics of Spin Glasses and In- 
formation Processing, Oxford University Press (2001). 

3 F. Merz and J. T. Chalker, Phys. Rev. B, 65, 054425 
(2002). 

4 P. W. Kasteleyn, J. Math. Phys., 4, 287 (1963). 

5 M. E. Fisher, 7, 10 (1966). 

6 L. Saul and M. Kardar, Phys. Rev. E, 48, R3221 (1993). 

7 L. Saul and M. Kardar, Nuc. Phys. B, 432, 641 (1994). 

8 A. Galluccio, M. Loebl, and J. Vondrak, Phys. Rev. Lett., 
84, 5924 (2000). 

9 R. M. F. Houtappel, Physica, 16, 391 (1950). 

I. Syozi, in Phase Transitions and Critical Phenomena, 
Domb and Green (eds), vol. 1 (1972). 



C. J. Colbourn, J. S. Provan, and D. Vertigan, Disc. Appl. 
Math., 60, 119. 

T. A. Feo and J. S. Provan, Oper. Res., 41, 572. 

D. J. Frank and C. J. Lobb, Phys. Rev. B, 37, 302 (1988). 
D. Archdeacon, C. J. Colbourn, I. Gitler, and J. S. Provan, 
Journal of Graph Theory, 33, 83 (2000). 

B. Derrida and L. de Seze, J. Physique, 43, 475 (1982). 
B. Derrida and J. Vannimenus, J. Phys. A, 15, L557 
(1982). 

A. George, SI AM J. Numer. Anal., 10, 345. 

R. J. Lipton, D. Rose, and R. E. Tarjan, SIAM J. Numer. 

Anal., 16, 346. 

I. Morgenstern and K. Binder, Phys. Rev. Lett., 43, 1615 
(1979). 

N. Robertson and P. D. Seymour, Journal of Combinato- 



12 



Type of model 


Distinct equalities 


Y parameters 


V parameters 


Ising 


4 


4 {F,K 01 ,K 02 ,K 03 ) 
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5 


4 
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I 3 


3q 2 - 2 


3q 2 - 2 



TABLE I: Degrees of freedom of Y and V sub-networks. For all cases except the standard Ising model, there are more equations 
than V parameters, indicating that the effective action contains 3-spin interactions and that exact Y-V transformations do not 
exist. 



Type of model 


Tystar-mcsh 
^ ' max 


Tyrmcsh-star 
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resistor networks, 
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OO 
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Ising model 


3 
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Ising with field 


2 


3 
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network reliability 


2 


2 
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transformations for various models. 
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